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ABSTRACT 

One-dimensional water-hammer modeling involves the solution of two coupled non-linear hyperbolic partial 
differential equations (PDEs). These equations result from applying the principles of conservation of mass and 
momentum to flow through a pipe, and usually the assumption that the speed at which pressure waves propagate 
through the pipe is constant. In order to solve these equations for the interested quantities (i.e. pressures and flow 
rates), they must first be converted to a system of ordinary differential equations (ODEs) by either approximating the 
spatial derivative terms with numerical techniques or using the Method of Characteristics (MOC). The MOC 
approach is ideal in that no numerical approximation errors are introduced in converting the original system of PDEs 
into an equivalent system of ODEs. Unfortunately this resulting system of ODEs is bound by a time step constraint 
so that when integrating the equations the solution can only be obtained at fixed time intervals. If the fluid system to 
be modeled also contains dynamic components (i.e. components that are best modeled by a system of ODEs), it 
may be necessary to take extremely small time steps during certain points of the model simulation in order to achieve 
stability and/or accuracy in the solution. Coupled together, the fixed time step constraint invoked by the MOC, and 
the occasional need for extremely small time steps in order to obtain stability and/or accuracy, can greatly increase 
simulation run times. As one solution to this problem, a method for combining variable step integration (VSI) 
algorithms with the MOC was developed for modeling water-hammer in systems with highly dynamic components. A 
case study is presented in which reverse flow through a dual-flapper check valve introduces a water-hammer event. 
The predicted pressure responses upstream of the check-valve are compared with test data. 

INTRODUCTION 

One dimensional incompressible transient flow problems fall into two categories: 1) bulk transient flow, and 
2) propagative transient flow. Historically, the propagative transient flow problem has been substantially more difficult 
to solve. This difficulty arises in the solution of the mathematical equations that model the differing assumptions in 
the previous two cases. This difference in assumptions boils down to the speed at which changes in pressure or 
pressure disturbances are propagated through the system, i.e. the effective speed of sound of the fluid in the system. 
For bulk transient flow, the speed of sound is assumed to be infinite; therefore an induced pressure disturbance at 
one point in the fluid system is felt everywhere in the system instantaneously. Conversely, for propagative transient 
flow the effective speed of sound of the fluid in the system is considered to be finite; therefore an induced pressure 
disturbance takes some finite amount of time to be felt everywhere in the system. In this case, the pressure 
disturbance propagates at the effective speed of sound of the fluid in the system. 

Although at first glance the difference in the two cases may seem somewhat trivial, mathematically the 
equations that result from these two assumptions are very different. At its root, modeling one dimensional 
incompressible transient flow is accomplished by solving the equations of conservation of mass and conservation of 
momentum for all necessary components that make up the fluid system. As it turns out, for bulk flow problems these 
equations are coupled sets of ODEs while for propagative flow they are coupled sets of hyperbolic PDEs. As luck 
would have it, at least at present, coupled sets of hyperbolic PDEs are some of the most difficult types of equations to 
solve. In fact, in order to solve these PDEs, it is necessary to first convert them to a system of ODEs by 
approximating all of the spatial derivatives by finite difference or finite element approximations, or by using the 
Method of Characteristics (MOC). 

In either the finite difference or finite element method one usually has to resort to an implicit scheme (i.e. a 
scheme which results in a system of coupled non-linear ODEs or algebraic equations that must be solved 
simultaneously) in order to obtain stability and/or accuracy in the solution. Solving large sets of non-linear equations, 
even for the one-dimensional problem, can result in a substantial increase in computational time. Particularly when 
the fluid system of interest contains other dynamic component models, i.e. models that occasionally must be 
integrated using very small time steps in order to preserve stability and/or accuracy in the solution. 

The MOC offers an attractive alternative to the finite difference and finite element methods. In solving the 
original PDEs, the MOC combined with Euler integration results in a system of algebraic equations in which all 
unknowns can be solved for explicitly (i.e. one at a time) given appropriate initial conditions. However, the MOC is 
bounded by a constraint that forces a fixed time step throughout the integration process. Typically for relatively 
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simple flow systems, the fixed integration time step resulting from the MOC solution scheme is not necessarily a 
disadvantage. The reason for this is the critical time step for the integration of the ODEs in the model is usually the 
integration time step dictated by using the MOC solution scheme. In other words, the rate at which the pressure 
waves propagate from one point to the next in the system is much faster than the rate of any other variables that are 
modeled. So, if your integration time step is small enough to capture pressure wave propagation, it is small enough 
to accurately capture everything else. However for more complex systems, i.e. systems in which component models 
consist of non-iinear ODEs, the integration time step might need to be very small over the course of a few integration 
steps in order to maintain the needed level of stability and/or accuracy in the solution. This very small time now 
becomes the critical time step. This time step might not be known in advance, and in this case would have to be 
determined by trial and error. Regardless, in these cases, a fixed time step integration routine could result in 
extremely long run times. The method developed in this paper couples a robust variable step integration (VSI) 
routine with the MOC in order to efficiently run models that do not produce numerical instabilities and oscillations in 
the solution. The paper also presents a comparison of model results with test data from a water flow experiment. 

RESULTS AND DISCUSSIONS 


MOC OVERVIEW FOR A SINGLE PIPE 


To ease explanation, the MOC solution applied to the water-hammer equations (neglecting gravity) is shown 

below. 
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Equations (1 ) and (2) are the original PDEs that represent one-dimensional water-hammer flow, while (3) and (4) 
constitute the set of ODEs that is produced when the MOC is applied to (1) and (2). (3b) and (4b) are the constraint 
or characteristic equations which link the pressure wave propagation speed (i.e. the effective speed of sound in the 
system) to the spatial grid size and integration time step. In solving (3a) we use upstream information constrained by 
(3b), and in solving (4a) we use downstream information constrained by (4b). It is clear from (3b) and (4b) that once 
the speed of sound and the grid size are chosen the integration time step becomes fixed. After these parameters are 
determined, (3a) and (4a) can be integrated (Euler integration is used in this case) subject to constraints (3b) and 
(4b); where each unknown is represented at some /* node along the one-dimensional space. 



C>M *= 0 


pa 2D ' ' 

(V,,.. - v ,D—(p u .» > = 0 


(5a) 

(6a) 

(5b) 

(6b) 


Above, (5b) and (6b) are the integrated form of (5a) and (6a) with the assumption that the third term in (5a) and (6a) 
will be held constant across each time step. At this point it is useful to rewrite (5b) and (6b) in a form that groups all 
known terms together, i.e. all terms at the current time. 
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Finally, we can combine (5c) and (6c) and rewrite them explicitly in terms of the unknowns, i.e. velocity and pressure. 
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Equations (7) and (8) allow for the solution of the unknowns at all internal nodes given initial conditions that satisfy 
(5b) and (6b). However, this does not hold at the boundaries (i.e. at the pipe inlet and exit) because (5d) has no 
meaning at the inlet and (6d) has no meaning at the exit. Therefore, in order to solve for the unknowns at the pipe 
inlet, (6c) must be solved in conjunction with its appropriate boundary condition relationships; while, at the pipe exit, 
(5c) must be solved in conjunction with its appropriate boundary condition relationships. 


VARIABLE STEP INTEGRATION AT THE BOUNDARIES 


Without upsetting the simplicity and efficiency of (7) and (8) for the internal nodes, our goal is to give validity 
to (5c) and (6c) at the boundaries even when the chosen integration time step does not abide by the constraint 
equations (3b) and (4b). Once accomplished, this allows the use of very robust pre-packaged ODE integrators to 
soive for variables of interest in complex fluid system models while using the MOC to pick up highly transient 
phenomenon in fluid lines. As it turns out, a method that works well is to linearly interpolate between nodal values of 
(5d) or (6d) (i.e. the C + or C' values) to a position that satisfies (3b) or (4b) for the time step chosen by the VSI. It is 
important to note that this applies at the boundaries only. In fact, variables at all internal nodes are only solved for 
when the integrated time matches a MOC time. Again, this falls back to (3b) and (4b). By fixing the grid size, (7) and 
(8) are only valid at MOC times. The following equations can be used to determine the distance from each boundary 
and therefore the spatial position where C + and C‘ are required so that (5c) and (6c) become valid at the boundaries 
at some time other than a MOC time. 


Ax 


i=F irst 


= +aA/„ 


Ax!:^' - -aA/„ 


(9a, 9b) 


Here, the change in x in (9a) is measured from the first node, while the change in x in (9b) is measured from the last 
node. It is important to note that all time step changes are measured from the last known MOC time. Once these 
spatial positions are calculated, the values for C + and C' can be determined by linearly interpolating from the two 
known values of C + or C' which bracket the position of interest. Using these interpolated values, (5c) and (6c) are 
valid at their appropriate boundaries. Therefore, with the appropriate boundary conditions the solution of all 
unknowns at the boundaries can be calculated for some time other than that dictated by the MOC. Of course there 
are a few practical considerations that must be taken into account when using these interpolated values to solve for 
variables at the pipe boundaries. The first is that the maximum time step allowed by the VSI should be limited to less 
than or equal to two times the change in time dictated by the MOC. This will insure that errors due to interpolation 
are held to a minimum. Also the VSI must integrate to all MOC times. This has to be done; it is only at the MOC 
times that the solution of the variables at the internal nodes may be obtained. Therefore the solution proceeds as 
follows: For some set of initial conditions begin to integrate. Use interpolated values of C + and C‘ to solve for all 
unknowns at the pipe boundaries as long as the VSI time step is less than MOC time step, again remembering that 
all time steps are measured from the last known MOC time. As soon as the VSI time is equal to MOC time, solve for 
all unknowns at the pipe boundaries and for all unknowns at the internal nodes within the pipe. Repeat the process 
until the desired solution time is reached. 
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Figure 1 

As an example consider the pipe inlet shown in Figure 1 . For any given time, C is known for all nodes, and 
the desire is to calculate all unknowns in the model at the next MOC time. Flowever, for stability and/or accuracy at 
the pipe boundaries the VSI takes a time step smaller than the MOC time step. In order for (6c) to be valid at this 
smaller time step a value for C' must be calculated that satisfies the constraint equation (4b) at the time chosen by 
the VSI. To achieve this, all that has to be done is to linearly interpolate between the C's at nodes 1 and 2 to the 
appropriate spatial position determined by (9a). If we use this new interpolated value of C in equation (6c), we have 
effectively shifted from the negative constraint at node 2, to the one labeled (A), and (6c) is valid for this time step. If 
VSI time step is larger than the MOC time step, the above procedure applies, but now we have to interpolate between 
C‘‘s at nodes 2 and 3 which causes a shift from the negative constraint at node 2 to the one labeled (B). Once the 
integration time reaches the first new MOC time, the time axis can be shifted down, and the procedure can be 
repeated. 

CASE STUDY : WATER FLOW CROSS-FEED EXPERIMENT 

The above method was used to model the test setup shown in Figure 2. During the initial phase of each 
test, water flowed from the low pressure tank through exit 1 , while also flowing through the cross-feed section and 
exit 2. Once a steady state condition was reached, the isolation vaive downstream of the high pressure tank was 
opened. This forced the check-valve closed so that during this phase water flowed from the low pressure tank 
through exit 1 and from the high pressure tank through exit 2. Over the course of testing, the following parameters 
were varied: 1) the difference in pressure between the low and high pressure tanks, 2) the steady state flow rates 
through the exits, and 3) the rate at which the isolation valve was opened. In this section, model results will be 
compared to test data for the pressure upstream of the check-valve for 2 of those tests. The two tests were selected 
because both produced water-hammer caused by flow reversal through the check-valve. The check-valve had an 
internal diameter of approximately 4 inches, and contained 2 spring loaded flaps connected at the center line. The 
flaps opened or closed depending on the pressure drop across the check-valve (see figure 3). 
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During the model build up phase it was determined that the check-valve model was critical in obtaining 
accurate predictions. This was because the pressure response of the system was very sensitive to the rate and time 
at which the check-valve closed. The two possible scenarios were: 1) the check-valve would close just as the flow 
rate in the cross-feed line went to zero, resulting in no dynamic pressure response; or 2) the check-valve would close 
after the flow rate reversed in the cross-feed line, possibly resulting in water-hammer. 


In order to incorporate a reasonable amount of physics the following check-valve model was used. 
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(10a) represents conservation of momentum for one check-valve flapper (see figure 3). The first term on the right 
hand side is the momentum exerted on the flapper due to the flowing fluid, while the second term is the momentum 
exerted by the flapper springs. Fp is a multiplier used to balance the steady state version of (10a) with manufacturer 
component level test data at the full open position. Rcp is the distance from the hinge to the center of pressure for 
each flapper. Ks is the spring constant and #ois the preload angle on the springs when the check-valve is closed. 
(10b) is conservation of momentum for the fluid flowing through the check valve. Aflap in (10a) and Kloss in (1 Ob) 
are both nonlinear functions of 9 . The above check-valve model coupled with standard models for tanks, valves, 
fittings, and an unaltered MOC line model (i.e. not coupled with a VSI) for all pipes in the system was initially used to 
simulate the water flow tests. 


Initially, during model simulation, a simple Euler fixed time step integration scheme was used. This decision 
was based on the following: 1) the normal MOC solution scheme needed a fixed integration time step, and 2) the 
choice of the grid size coupled with the constraint equations forced the time step to be on the order of 1 E-4 seconds. 
However, in coupling the check-valve physics with the rest of the model, the model became unstable. The Euler 
integration scheme did not work, even with time steps as small as 1 E-6 seconds. The problems appeared to occur 
just as the check valve was starting to open and when it was going from partially opened to closed. An integration 
algorithm and time step sensitivity study was performed on the check-valve model for various sinusoidal pressure 
drop inputs. As a ground rule, the integration algorithms were restricted to the Euler and Gear 


I 



Spring Force 
Moment 


i 

Figure 3 


schemes. Fluid system models are known to be numerically stiff; therefore the Runge-Kutia scheme was not 
evaluated as an option. The Gear algorithm was a VSI, and was used to determine at what points of check-valve 
operation (i.e. going from closed to open, or opened to closed) the time step had to be dramatically reduced in order 
to obtain stability and/or accuracy in the solution. Using the Euler integration scheme the time step had to be driven 
down to approximately IE-8 seconds in order to achieve stability and/or accuracy as measured by the Gear 
algorithm. In other words, because the Gear algorithm was a VSI, it was assumed that it maintained a certain level of 
accuracy in its solution of the check-valve model as dictated by its user defined tolerance settings. Therefore the 
Euler solution approximated the Gear solution when a time step of 1 E-8 seconds was used. This fixed time step was 
not feasible; model run times would have been extremely long. The same sinusoidal pressure drop input runs were 
made using Gear’s algorithm. This method worked well and was efficient. However, even this algorithm was forced 
to time steps on the order or IE-6 seconds at some tough points during the integration process. These results led to 
the conclusion that a VSI would have to be coupled with the MOC solution scheme in order to achieve stable and 
accurate results in a relatively timely manner. 

MODEL VERSUS TEST DATA 


Figures (4a) and (5a) show the pressure just upstream of the check-valve during the transient response 
caused by the isolation valve opening and the check-valve closing. The test parameters for the cases shown in 
figures (4) and (5) are identical, except for the rate at which the isolation valve opens. Both cases have a steady 
state flow rate through the cross-feed line of 500 gpm with a difference between the high and low tank pressures of 
1 5 psig. The time it takes the isolation valve to open for the figure (4) case is 3.0 sec, while its 0.5 sec for the figure 
(5) case. By inspection of both figures (4a) and (5a), it can be seen that the method presented herein predicts the 
initial pressure surge caused by the isolation valve cracking open very closely. The pressure at times between the 
isolation valve opening and the check-valve closing is also predicted well. As the check valve closes, the pressure 
initially decreases as shown in figures (4) and (5). The decrease is caused by the flow direction in the cross-feed line 
as the check-valve goes closed. In both cases the flow had reversed before the check-valve closed. Figures (4b) 
and (5b) show an expanded view of the pressure response around the time of check-valve closure. The method also 
predicted this pressure response closely for the first few cycles. There is a lack of attenuation in the model results 
after the first few pressure cycles. This could be caused by energy absorption in the physical system that is not 
represented in the model physics as presented here, and/or by a dependency of frictional loss on the rate of change 
of velocity through both time and space, i.e. frequency dependent friction. It should be noted that model results 
presented here used one form of a frequency dependent friction term. The pressure signature after check-valve 
closure produced without the frequency dependent friction term took much longer to attenuate than what is presented 
here. 
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Figure 5a 
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Figure 5b 


NOMENCLATURE 


Afiap 

ft 2 

area of check-valve flapper 

a 

ft/sec 

effective speed of sound 

D 

ft 

diameter 

f p 

none 

moment muitipiier 

f 

none 

friction factor 

9c 

lbm-ft/lbf-sec 2 

conversion factor 

Kioss 

none 

check-valve loss coefficient 

Ks 

Ibf-ft/rad 

check-valve spring constant 

1 

lbf-ft-sec 2 

mass moment of inertia for check-valve flapper 

p 

Ibf/ft 2 

pressure 

Pip 

Ibf/ft 2 

pressure at check-valve inlet 

P out 

Ibf/ft 2 

pressure at check-valve exit 

Rco 

ft 

distance from check-valve hinge to flapper center of pressure 

t 

sec 

time 

V 

ft/sec 

velocity 

X 

ft 

spatial position 

p 

Ibm/ft 2 

density 

e 

rad 

check-valve flapper angle 

00 

rad 

check-valve flapper pre-load angle 


SUMMARY AND CONCLUSIONS 

A method is presented that allows the modeling of fluid systems when propagative transient flow is important 
and the system also contains highly dynamic components. The method couples the MOC solution of the one- 
dimensional water-hammer equations with a VSI routine. This allows for the stability, accuracy, and speed benefits 
gained by using the MOC to solve these equations. It also affords the ability to dramatically reduce the integration 
time step in order to maintain stability and/or accuracy in the solution of other highly dynamic components in the fluid 
system. A case study is presented which compares test data to model results for a given pressure location. The fluid 
system causes transient pressure responses by an opening isolation valve that forces flow reversal in a fluid line 
which in turn forces a 4 inch dual flapper check valve closed. The method presented herein predicted these transient 
responses reasonably accurately. 
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